package defectsim;

import defectsim.DefectFactory.DefectType;
import repast.simphony.random.RandomHelper;
import repast.simphony.space.continuous.NdPoint;

public class EdgeDislocation extends Dislocation {

	private Material mat;
	private boolean shortRangeInteraction;
	private boolean dipoleFormation;
	private double passingShearStress;
	private Iterable<Dislocation> longRangeNeighbours;
	private Iterable<Dislocation> midRangeNeighbours;
	private Iterable<Dislocation> shortRangeNeighbours;
	private double glideStress;
	private double longRangeStress;
	private double effectiveStress;
	private double orowanStress;
	private double acceleration;
	double longPerceptionRadius;
	double shortPerceptionRadius;
	double midPerceptionRadius;
	double[] stressForRK4 = new double[4];
	
	SimEnvironment env = SimEnvironment.getInstance();

	@Override
	public void percept() {
		mat = env.getCurrentMaterial();
		// Stochastic cutoff radius
		longPerceptionRadius = RandomHelper.nextDoubleFromTo(env.getMinCutOffRadius(),
				env.getMaxCutOffRadius());
		// Deterministic cutoff radius
//		longPerceptionRadius = env.getMaxCutOffRadius();
		shortPerceptionRadius = env.getAnnihilationLength() * 1e-3; // in microns
		midPerceptionRadius = RandomHelper.nextDoubleFromTo(0.5 * env.getDipoleLength(),
				env.getDipoleLength()) * 1e-3; // in microns
		longRangeNeighbours = env.getAllNeighbours(this);
		midRangeNeighbours = env.getNeighbours(this, midPerceptionRadius);
		shortRangeNeighbours = env.getNeighbours(this, shortPerceptionRadius);
		// Revision for RK4 integration
//		longRangeStress = Util.calculateTotalLongRangeStress(this, longRangeNeighbours, env, longPerceptionRadius);
		stressForRK4 = Util.calculateTotalLongRangeStressForRK4(this, longRangeNeighbours, env, longPerceptionRadius);
		longRangeStress = stressForRK4[0] - env.getCurrentExternalStress() + env.getCurrentFrictionStress();
		effectiveStress = env.getCurrentExternalStress() + longRangeStress;
		glideStress = effectiveStress - env.getCurrentFrictionStress();
//		double mys = calculateStress();
//		System.out.printf("id: %10s, LR: %f, Ef: %f, Gs: %f\n", this.toString(), longRangeStress, effectiveStress, glideStress);
//		System.out.printf("LR: %f, Mys: %f\n", longRangeStress, mys);
		this.setStress(effectiveStress);
		passingShearStress = mat.getShearModulus() * this.getUnitBurgersVector() * mat.getBurgersVector() /
		(8 * Math.PI * (1 - mat.getPoissonRatio()) * mat.getSlipPlaneSpacing());

		// Orowan stress criterion Eq. 27 from "Dislocation Dynamics I. A Proposed...", R.J. Amodeo & N.M.Ghoniem, 1990
		orowanStress = Math.abs(2 * mat.getShearModulus() * this.getUnitBurgersVector() * mat.getBurgersVector() / mat.getInterobstacleSpacing());
	}

	@Override
	public void mentalUpdate(){
		// If there is at least one dislocation in the vicinity of this dislocation 
		// it encounters short range interaction, e.g. annihilation or junction formation
		if (shortRangeNeighbours != null) shortRangeInteraction = true;

		// If there is at least one dislocation in the mid range of this dislocation 
		// there is a possibility for dipole formation. The second criteria for dipole formation
		// is that the effective stress on the dislocation should be less than passing shear stress
		if (midRangeNeighbours != null && (effectiveStress < passingShearStress)) dipoleFormation = true;

		// Mobility rule
		double dt = env.getTimeStep();
		double v = this.getVelocity();
		double d = 0;
		if (Math.abs(effectiveStress) > orowanStress) {
//			System.out.println("Orowan stress criterion is satisfied."+ effectiveStress + " " + orowanStress);
			env.setMultiplication(true);
		}
		if ((Math.abs(effectiveStress) > env.getCurrentFrictionStress()) && !this.isDipole()) {
			this.setMobility(true); // mobile
			double force = stressForRK4[0] * this.getUnitBurgersVector() * mat.getBurgersVector();
			// Velocity calculation method 1
//			// Velocity Verlet algorithm is applied for time integration of equations of motion
//			d = v * dt + acceleration * dt * dt / 2.0; 
//			double vtemp = v + acceleration * dt / 2.0;
//			acceleration = force * 3.38e7;
//			v = vtemp + acceleration * dt / 2.0;
//			if(v == Double.POSITIVE_INFINITY){
//				System.out.println("infinity found!");
//				System.out.println(this.toString() + "Force: " + force);
//				System.out.println(this.toString() + "Acc.: " + acceleration);
//				System.out.println(this.toString() + "Dist: " + d);
//				System.out.println(this.toString() + "Vtemp: " + vtemp);
//			}
			// Euler integration
			
			// Velocity calculation method 2 (The first that I tried)
//			double nanometerTOmeter = 1e-09;
//			double v0 = 1e-2;
//			v = (v0 / 2700) * (force / (this.getUnitBurgersVector() * 0.249 * nanometerTOmeter));
			
			// Velocity calculation method 3
			// Velocity in meter/second
			v = force / mat.getDampingCoefficient();
			
			// Runge Kutta 4th order integration
			// Velocity calculation method 4 based on RK4 integration
//			double r = (1 / mat.getDampingCoefficient()) * this.getUnitBurgersVector() * mat.getBurgersVector() * (dt / 6) * 
//				(stressForRK4[0] + 2 * stressForRK4[1] + 2 * stressForRK4[2] + stressForRK4[3]) ;
//			v = r / dt;
			
			// Distance in micrometer
			// The formula below should be used for method 1, 2 & 3
			d = env.getTimeStep() * v * 1e6;
		}
		else {
			this.setMobility(false); // immobile
			v = 0;
		}
		this.setVelocity(v);
		this.setDistance(d);
	}

	@Override
	public void act() {
		boolean annihilated = false;
		if (shortRangeInteraction){
			// check for annihilation
			for (Dislocation d : shortRangeNeighbours) {
				double dir = d.getDirection();
				double uBurgersVec = d.getUnitBurgersVector();
				if (dir == this.getDirection()) 
					if (uBurgersVec == -1.0 * this.getUnitBurgersVector()){
						env.annihilate(d);
						env.annihilate(this);
						annihilated = true;
					}
					else{
						this.setUnitBurgersVector(this.getUnitBurgersVector() + d.getUnitBurgersVector());
						env.annihilate(d);
					}

			}
		}
		if (dipoleFormation && !this.isDipole() && !annihilated){
			// check for dipole formation
			NdPoint p1 = env.getLocation(this);
			for (Dislocation d : midRangeNeighbours) {
				if(!d.isDipole()){
					NdPoint p2 = env.getLocation(d);
					// check whether the dislocation is annihilated before 
					if ((p2 != null) && (d.getUnitBurgersVector() == -1.0 * this.getUnitBurgersVector())){
						double dist = Util.calculateDistance(p1, p2) * 1e3; // in nm
						double teta = Util.calculateAngle(p1, p2) * 180 / Math.PI; // in degrees
						if ((dist > env.getAnnihilationLength()) && 
								((teta > 30.0 && teta < 60.0) ||
								 (teta > 120.0 && teta < 150.0) ||
								 (teta > 210.0 && teta < 240.0) ||
								 (teta > 280.0 && teta < 330.0))){
							this.setDipole(true);
							d.setDipole(true);
							this.setMobility(false);
							this.setVelocity(0.0);
							this.setDipoleMate(d);
							d.setDipoleMate(this);
							break;
						}
					}
				}
			}
		}
		if (this.isDipole() && !annihilated){
			// destruct dipole formation if effective stress is greater than passing shear stress
			if (effectiveStress > passingShearStress){
				this.setDipole(false);
				this.setMobility(true);
				if (this.getDipoleMate() != null)
					this.getDipoleMate().setDipole(false);
			}
		}
		if (this.isMobile() && !annihilated)
			env.moveByVector(this, this.getDistance(), Math.toRadians(this.getDirection()));


	}

	private double calculateStress(){
		double stress = 0.0;
		double length = env.getSpaceLength();
		NdPoint p1 = env.getLocation(this);
		for (Dislocation dj : longRangeNeighbours) {
			NdPoint p2 = env.getLocation(dj);
			double xx = p2.getX() - p1.getX();
			double yy = p2.getY() - p1.getY();
			double dx = xx - length * Math.round(xx / length);
			double dy = yy - length * Math.round(yy / length);
			double r = Math.sqrt(dx * dx + dy * dy);
			double delta = StrictMath.atan2(dy, dx);
			while(0 > delta) {
				delta += 2.0 * Math.PI;
			}
			double teta = delta - Math.toRadians(this.getDirection());
			double gamma = delta - Math.toRadians(dj.getDirection());
			if (r > 0 && r <= longPerceptionRadius){
				r *= 1000;
				
				double A = (mat.getShearModulus() * dj.getUnitBurgersVector() * mat.getBurgersVector()) / 
				(2 * Math.PI * (1 - mat.getPoissonRatio()));

				stress = stress + A * Math.cos(teta) * Math.cos(2 * gamma)/ r;
			}
		}
		return stress;
	}
	
	
	public EdgeDislocation(){
		super.setDefectType(DefectType.EdgeDislocation);
		super.setMobility(true); // mobile
		super.setDipole(false); // monopolar
		super.setVelocity(0.0);
		dipoleFormation = false; // no possible dipole formation
		shortRangeInteraction = false; // no short range interaction
	}

	
	public int isEdgeDislocation(){
		return 1;
	}
}
